##Peacekeepers without helmets: How violence shapes local peacebuilding by civilian peacekeepers
##Hannah Smidt & Allard Duursma
## Replication for Appendix I

library(rgdal)
library(dplyr)
library(splm)
library(rgeos)
library(sf)
library(conleyreg)

# Set working directory
setwd("SET YOUR OWN WORKING DIRECTORY")

# Read final data
data <- read.csv("2022_11_07_dataFinal_max_agg_MonEv.csv")


# Map of CAR
folder <- "./caf_admbnda_adm3_200k_sigcaf_reach_itos_v2/"
map      <- readOGR(dsn = folder, layer = "caf_admbnda_adm3_200k_sigcaf_reach_itos_v2")
map_poly <- sp::SpatialPolygons(map@polygons)

# Make spatial weights matrix
nb <- spdep::poly2nb(map_poly)
w <-  spdep::nb2mat(nb, style="B", zero.policy=TRUE)
wlist <- spdep::nb2listw(nb, style="B", zero.policy=TRUE)


# Order data by time period
data$month <- as.numeric(data$month)
data$year <- as.numeric(data$year)
data2 <- data %>% arrange(year,month,admin1Name,admin2Name,admin3Name,admin3RefN)
data2$id <- apply(data2[c("admin1Name", "admin2Name", "admin3Name", "admin3RefN")],1,function(x) paste(x, collapse="") )
data2$id <- as.numeric(as.factor(data2$id))
data2$day<-1
data2$yearmo <- as.Date(apply(data2[c("year", "month", "day")],1,function(x) paste(x, collapse="/") ), "%Y/%m/%d" )
data2$yearmo <- zoo::as.yearmon(data2$yearmo)

# Get rid of observations outside time frame
data2 <- subset(data2, !(year==2018 & month>=11))

prop.table(table(data2$LocalPeacebuild_AnyAssistance))
sd(data2$LocalPeacebuild_AnyAssistance)

# Make panel and check balance
pdata_splm <- plm::pdata.frame( data2, index = c("id", "yearmo") )
plm::is.pbalanced(pdata_splm)




#################################################################################
## Conley Standard errors accounting for a) serial correlation ##################
## and spatial dependence #######################################################
#################################################################################


# Make simple features data frame
centroids <- SpatialPointsDataFrame(gCentroid(map, byid = TRUE, id = map@data$OBJECTID_1),  map@data ) 

data3 <- merge(data2, centroids, by=c("admin1Name", "admin2Name", "admin3Name", "admin3RefN"), all.y=T)

coordinates(data3) <- ~ x + y
data3 <- sf::st_as_sf(data3)
st_crs(data3) = 4326


m1a <- conleyreg(LocalPeacebuild_AnyAssistance ~ ACLED_viol_any, data = data3, dist_cutoff = 1000)


m1b <- conleyreg(LocalPeacebuild_AnyAssistance ~ ACLED_viol_any + roadDensity + distToCapital +
                   Shape_Area_adm2 + Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize +
                   ACLED_viol_any_3m_MA + anyCheckpoint, data = data3, dist_cutoff = 1000)


m1c <- conleyreg(LocalPeacebuild_AnyAssistance ~ ACLED_viol_any + FoodInsecurity +
                   ACLED_viol_any_3m_MA + anyCheckpoint | id, data = data3, dist_cutoff = 1000)

m1a
m1b
m1c

 



######################################################
## Other Spatial error models (NOT IN APPENDIX) ###
######################################################


m1a <- spml(formula = (LocalPeacebuild_AnyAssistance ~ ACLED_viol_any)
            , data = pdata_splm, listw = wlist, model = "pooling"
            , cl=pdata_splm$id, spatial.error=c("b") )
summary(m1a)

m1b <- spml(formula = (LocalPeacebuild_AnyAssistance ~ ACLED_viol_any + roadDensity + distToCapital +
                         Shape_Area_adm2 + Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize +
                         ACLED_viol_any_3m_MA + anyCheckpoint)
            , data = pdata_splm, listw = wlist, model = "pooling"
            , spatial.error=c("b") )
summary(m1b)

m1c <- spml(formula = (LocalPeacebuild_AnyAssistance ~ ACLED_viol_any + FoodInsecurity +
                         ACLED_viol_any_3m_MA + anyCheckpoint)
            , data = pdata_splm, listw = wlist, model = "within"
            , spatial.error=c("b") )
summary(m1c)

